Note: The description became a little longer than expected. Do you know a readable implementation of this algorithm using this mesh? Please let me know!

I'm trying to implement Catmull-Clark subdivision using Matlab (because later on the results have to be compared with some other stuff already implemented in Matlab).

First try was with a Vertex-Face mesh, the algorithm works but it is of course not very efficient (since you need neighbouring information for edges and faces). Therefore, I'm now using a Half-Edge mesh (info), see also the paper of Lutz Kettner. Wikipedia link to the idea behind Catmull-Clark SDV: Wiki.

My problem lies in finding the Twin HalfEdges, I'm just not sure how to do this. Below I'm describing my thoughts on the implementation, trying to keep it concise.

Half-Edge mesh (using indices to Vertices/HalfEdges/Faces):

Vertex (x,y,z,Outgoing_HalfEdge)
HalfEdge (HeadVertex (or TailVertex, which one should I use), Next, Face, Twin).
Face (HalfEdge)

To keep it simple for now, assume that every face is a quadrilateral. The actual mesh is a list of Vertices, HalfEdges and Faces. The new mesh will consist of NewVertices, NewHalfEdges and NewFaces, like this (note: Number_... is the number of ...):

NumberNewVertices: Number_Faces + Number_HalfEdges/2 + Number_Vertices
NumberNewHalfEdges: 4 * 4 * NumberFaces
NumberNewfaces: 4 * NumberFaces

Catmull-Clark:

Find the FacePoint (centroid) of each Face:
--> Just average the x,y,z values of the vertices, save as a NewVertex. 
Find the EdgePoint of each HalfEdge:
--> To prevent duplicates (each HalfEdge has a Twin which would result in the same HalfEdge)
--> Only calculate EdgePoints of the HalfEdge which has the lowest index of the Pair.
Update old Vertices

Ok, now all the new Vertices are calculated (however, their Outgoing_HalfEdge is still unknown). Next step to save the new HalfEdges and Faces. This is the part causing me problems!

Loop through each old Face, there are 4 new Faces to be created
  (because of the quadrilateral assumption)
First create the 4 new HalfEdges per New Face, starting at the FacePoint to the Edgepoint
Next a new HalfEdge from the EdgePoint to an Updated Vertex
Another new one from the Updated Vertex to the next EdgePoint
Finally the fourth new HalfEdge from the EdgePoint back to the FacePoint.

The HeadVertex of each new HalfEdge is known, the Next HalfEdge too. The Face is also known (since it is the new face you're creating!). Only the Twin HalfEdge is unknown, how should I know this?

By the way, while looping through the Vertices of the new Face, assign the Outgoing_HalfEdge to the Vertices. This is probably the place to find out which HalfEdge is the Twin.

Finally, after the 4 new HalfEdges are created, save the Face with the HalfVertex index the last newly created HalfVertex.

I hope this is clear, if needed I can post my (obviously not-yet-finished) Matlab code.